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Abstract 

A general formulation was developed to represent material models for applications in dynamic 
loading. Numerical methods were devised to calculate response to shock and ramp compression, and 
ramp decompression, generalizing previous solutions for scalar equations of state. The numerical 
methods were found to be flexible and robust, and matched analytic results to a high accuracy. 
The basic ramp and shock solution methods were coupled to solve for composite deformation 
paths, such as shock-induced impacts, and shock interactions with a planar interface between 
different materials. These calculations capture much of the physics of typical material dynamics 
experiments, without requiring spatially-resolving simulations. Example calculations were made of 
loading histories in metals, illustrating the effects of plastic work on the temperatures induced in 
quasi-isentropic and shock-release experiments, and the effect of a phase transition. 
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I. INTRODUCTION 



The continuum representation of matter is widely used for material dynamics in sci- 
ence and engineering. Spatially-resolved continuum dynamics simulations are the most 
widespread and familiar, solving the initial value problem by discretizing the spatial domain 
and integrating the dynamical equations forward in time to predict the motion and defor- 
mation of components of the system. This type of simulation is used, for instance, to study 
hypervelocity impact problems such as the vulnerability of armor to projectiles p], the 
performance of satellite debris shields g and the impact of meteorite, with planets, notably 
the formation of the moon [4|]. The problem can be divided into the dynamical equations 
of the continuum, the state field of the components s(r), and the inherent properties of 
the materials. Given the local material state s, the material properties allow the stress r 
to be determined. Given the stress field r(r) and mass density field p(r), the dynamical 
equations describe the fields of acceleration, compression, and thermodynamic work done 
on the materials. 

The equations of continuum dynamics describe the behavior of a dynamically deforming 
system of arbitrary complexity. Particular, simpler deformation paths can be described more 
compactly by different sets of equations, and solved by different techniques than those used 
for continuum dynamics in general. Simpler deformation paths occur often in experiments 
designed to develop and calibrate models of material properties. These paths can be regarded 
as different ways of interrogating the material properties. The principal examples in material 
dynamics are shock and ramp compression 0,0]. Typical experiments are designed to induce 
such loading histories and measure or infer the properties of the material in these states 
before they are destroyed by release from the edges or by reflected waves. 

The development of the field of material dynamics was driven by applications in the 
physics of hypervelocity impact and high explosive systems, including nuclear weapons ?J ■ 
In the regimes of interest, typically components with dimensions ranging from millime- 
ters to meters and pressures from 1 GPa to 1 TPa, material behavior is dominated by the 
scalar equation of state (EOS): the relationship between pressure, compression (or mass 
density), and internal energy. Other components of stress (specifically shear stresses) are 
much smaller, and chemical explosives react promptly so can be treated by simple mod- 
els of complete detonation. EOS were developed as fits to experimental data, particularly 
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to series of shock states and to isothermal compression measurements [8|. It is relatively 
straightforward to construct shock and ramp compression states from an EOS algebraically 
or numerically depending on the EOS, and to fit an EOS to these measurements. More 
recently, applications and scientific interest have grown to include a wider range of pressures 
and time scales, such as laser-driven inertial confinement fusion jj], and experiments are 
designed to measure other aspects than the EOS, such as the kinetics of phase changes, con- 
stitutive behavior describing shear stresses, incomplete chemical reactions, and the effects of 
microstructure, including grain orientation and porosity. Theoretical techniques have also 
evolved to predict the EOS with ~1% accuracy [10|] and elastic contributions to shear stress 



with slightly poorer accuracy 



A general convention for representing material states is described, and numerical methods 
are reported for calculating shock and ramp compression states from general representations 
of material properties. 



II. CONCEPTUAL STRUCTURE FOR MATERIAL PROPERTIES 

The desired structure for the description of the material state and properties under dy- 
namic loading was developed to be as general as possible with respect to the types of material 
or models to be represented in the same framework, and designed to give the greatest amount 
of commonality between spatially-resolved simulations and calculations of shock and ramp 
compressions. 

In condensed matter on sub-microsecond time scales, heat conduction is often too slow to 
have a significant effect on the response of the material, and is ignored here. The equations 
of non-relativistic continuum dynamics are, in Lagrangian form, i.e. along characteristics 
moving with the local material velocity u(r), 

Dp(f,t) 



Dt 

Du(f, t) 1 



Dt p{f,t) 
De(f, t) 



—p(r,t)divu(r,t) (1) 
div r(f,t) (2) 



Dt 



r(r, t)gradu(r, t)\\ (3) 



where p is the mass density and e the specific internal energy. Changes in e can be related 
to changes in the temperature T through the heat capacity. The inherent properties of 
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each material in the problem are described by its constitutive relation or equation of state 
r(s). As well as experiencing compression and work from mechanical deformation, the local 
material state s(f,t) can evolve through internal processes such as plastic flow. In general, 



which can also include the equations for dp/dt and de/dt. Thus the material properties must 
describe at a minimum r(s) and s[s(f, t), U(r, t)] for each material. If they also describe T(s), 
the conductivity, and s(e), then heat conduction can be treated. Other functions may be 
needed for particular numerical methods in continuum dynamics, such as the need for wave 
speeds (e.g. the longitudinal sound speed), which are needed for time step control in explicit 
time integration. Internally, within the material properties models, it is desirable to re-use 
software as much as possible, and other functions of the state are therefore desirable to allow 
models to be constructed in a modular and hierarchical way. Arithmetic manipulations must 
be performed on the state during numerical integration, and these can be encoded neatly 
using operator overloading, so the operator of the appropriate type is invoked automatically 
without having to include 'if-then-else' structures for each operator as is the case in non- 
object-oriented programming languages such as Fortran- 77. For instance, if s is calculated 
in a forward-time numerical method then changes of state are calculated using numerical 
evolution equations such as 



Thus for a general state s and its time derivative s, which has an equivalent set of compo- 
nents, it is necessary to multiply a state by a real number and to add two states together. 
For a specific software implementation, other operations may be needed, for example to 
create, copy, or destroy a new instance of a state. 

The attraction of this approach is that, by choosing a reasonably general form for the 
constitutive relation and associated operations, it is possible to separate the continuum 
dynamics part of the problem from the inherent behavior of the material. The relations 
describing the properties of different types of material can be encapsulated in a library form 
where the continuum dynamics program need know nothing about the relations for any spe- 
cific type of material, and vice versa. The continuum dynamics programs and the material 
properties relations can be developed and maintained independently of each other, provided 
that the interface remains the same (Tabled]). This is an efficient way to make complicated 




U = grad u(f, t) 



(4) 



s(t + St) = s(t) + Sts. 



(5) 
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material models available for simulations of different types, including Lagrangian and Eule- 
rian hydrocodes operating on different numbers of dimensions, and calculations of specific 
loading or heating histories such as shock and ramp loading discussed below. Software in- 
terfaces have been developed in the past for scalar EOS with a single structure for the state 
fjj], but object-oriented techniques make it practical to extend the concept to much more 
complicated states, to combinations of models, and to alternative types of model selected 
when the program is run, without having to find a single super-set state encompassing all 
possible states as special cases. 

A very wide range of types of material behavior can be represented with this formalism. 
At the highest level, different types of behavior are characterized by different structures for 
the state s (Table HT|) . For each type of state, different specific models can be defined, such 
as perfect gas, polytropic and Griineisen EOS. For each specific model, different materials 
are represented by choosing different values for the parameters in the model, and different 
local material states are represented through different values for the components of s. In the 
jargon of object-oriented programming, the ability to define an object whose precise type 
is undetermined until the program is run is known as polymorphism. For our application, 
polymorphism is used at several levels in the hierarchy of objects, from the overall type of a 
material (such as 'one represented by a pressure-density-energy EOS' or 'one represented by 
a deviatoric stress model') through the type of relation used to describe the properties of that 
material type (such as perfect gas, polytropic, or Griineisen for a pressure-density-energy 
EOS, or Steinberg-Guinan 13] or Preston- Tonks- Wallace 14] for a deviatoric stress model), 
to the type of general mathematical function used to represent some of these relations (such 
as a polynomial or a tabular representation of j(p) in a polytropic EOS) (Table IITT|) . States 
or models may be defined by extending or combining other states or models - this can be 
implemented using the object-oriented programming concept of inheritance. Thus deviatoric 
stress models can be defined as an extension to any pressure-density-energy EOS (they are 
usually written assuming a specific type, such as Steinberg's cubic Griineisen form) , homo- 
geneous mixtures can be defined as combinations of any pressure-density-temperature EOS, 
and heterogeneous mixtures can be defined as combinations of materials each represented 
by any type of material model. 

Trial implementations have been made as libraries in the C++ and Java programming 
languages 15]. The external interface to the material properties was general at the level 



5 



of representing a generic material type and state. The type of state and model were then 
selected when programs using the material properties library were run. In C++, objects 
which were polymorphic at run time had to be represented as pointers, requiring additional 
software constructions to allocate and free up physical memory associated with each object. 
It was possible to include general re-usable functions as polymorphic objects when defining 
models: real functions of one real parameter could be polynomials, transcendentals, tabular 
with different interpolation schemes, piecewise definitions over different regions of the one 
dimensional line, sums, products, etc; again defined specifically at run time. Object-oriented 
polymorphism and inheritance were thus very powerful techniques for increasing software 
re-use, making the software more compact and more reliable through the greater use of 
functions which had already been tested. 

Given conceptual and software structures designed to represent general material proper- 
ties suitable for use in spatially-resolved continuum dynamics simulations, we now consider 
the use of these generic material models for calculating idealized loading paths. 

III. IDEALIZED ONE-DIMENSIONAL LOADING 

Experiments to investigate the response of materials to dynamic loading, and to calibrate 
parameters in models of their behavior, are usually designed to apply as simple a loading 
history as is consistent with the transient state of interest. The simplest canonical types of 
loading history are shock and ramp . Methods of solution are presented for calculating 
the result of shock and ramp loading for materials described by generalized material models 
discussed in the previous section. Such direct solution removes the need to use a time- 
and space-resolved continuum dynamics simulation, allowing states to be calculated with 
far greater efficiency and without the need to consider and make allowance for attributes of 
resolved simulations such as the finite numerical resolution and the effect of numerical and 
artificial viscosities. 

A. Ramp compression 

Ramp compression is taken here to mean compression or decompression. If the material 
is represented by an inviscid scalar EOS, i.e. ignoring dissipative processes and non-scalar 
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effects from elastic strain, ramp compression follows an isentrope. This is no longer true 
when dissipative processes such as plastic heating occur. The term 'quasi-isentropic' is 
sometimes used in this context, particularly for shockless compression; here we prefer to 
refer to the thermodynamic trajectories as adiabats since this is a more appropriate term: 
no heat is exchanged with the surroundings on the time scales of interest. 

For adiabatic compression, the state evolves according to the second law of thermody- 
namics, 

de = T dS — pdv (6) 
where T is the temperature and S the specific entropy. Thus 

e = TS — pv = TS , (7) 

or for a more general material whose stress tensor is more complicated than a scalar pressure, 

de = T dS + r n dv =>- e = TS H — " (8) 

P 

where r n is the component of stress normal to the direction of deformation. The velocity 
gradient was expressed through a compression factor 77 = p' j p and a strain rate e. In all 
ramp experiments used in the development and calibration of accurate material models, 
the strain has been applied uniaxially. More general strain paths, for instance isotropic or 
including a shear component, can be treated by the same formalism, and that the working 
rate is then a full inner product of the stress and strain tensors. 

The acceleration or deceleration of the material normal to the wave as it is compressed 
or expanded adiabatically is 

Du __ j dr n 

Dv~ V dv' 1 J 

from which it can be deduced that 

d-p = i (10) 

where q is the longitudinal wave speed. 

As with continuum dynamics, internal evolution of the material state can be calculated 
simultaneously with the continuum equations, or operator split and calculated periodically 



at constant compression 16j. The results are the same to second order in the compression 
increment. Operator-splitting allows calculations to be performed without an explicit en- 
tropy, if the continuum equations are integrated isentropically and dissipative processes are 
captured by internal evolution at constant compression. 



Operator-splitting is desirable when internal evolution can produce highly nonlinear 
changes, such as reaction from solid to gas: rapid changes in state and properties can 
make numerical schemes unstable. Operator-splitting is also desirable when the integration 
time step for internal evolution is much shorter than the continuum dynamics time step. 
Neither of these considerations is very important for ramp compression without spatial res- 
olution, but operator-splitting was used as an option in the ramp compression calculations 
for consistency with continuum dynamics simulations. 

The ramp compression equations were integrated using forward-time Runge-Kutta nu- 
merical schemes of second order. The fourth order scheme is a trivial extension. The 
sequence of operations to calculate an increment of ramp compression is as follows: 

1. Time increment: 

I In 77 1 . . 

St=-^—r 11 (11) 

2. Predictor: 

s(t + 8t/2)=s(t) + ^s m (s(t),e) (12) 

3. Corrector: 

s(t + St) = s(t) + 5ts m (s(t + 5t/2), e) (13) 

4. Internal evolution: 

/t+st 
Si(s(t'),e)dt' (14) 

where s m is the model-dependent state evolution from applied strain, and is internal 
evolution at constant compression. 

The independent variable for integration is specific volume v or mass density p; for 
numerical integration finite steps are taken in p and v. The step size Ap can be controlled so 
that the numerical error during integration remains within chosen limits. A tabular adiabat 
can be calculated by integrating over a range of v or p, but when simulating experimental 
scenarios the upper limit for integration is usually that one of the other thermodynamic 
quantities reaches a certain value, for example that the normal component of stress reaches 
zero, which is the case on release from a high pressure state at a free surface. Specific 
end conditions were found by monitoring the quantity of interest until bracketed by a finite 
integration step, then bisecting until the stop condition was satisfied to a chosen accuracy. 
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During bisection, each trial calculation was performed as an integration from the first side 
of the bracket by the trial compression. 



B. Shock compression 

Shock compression is the solution of a Riemann problem for the dynamics of a jump 
in compression moving with constant speed and with a constant thickness. The Rankine- 
Hugoniot (RH) equations [5J describing the shock compression of matter are derived in 
the continuum approximation, where the shock is a formal discontinuity in the continuum 
fields. In reality, matter is composed of atoms, and shocks have a finite width governed by 
the kinetics of dissipative processes - at a fundamental level, matter does not distinguish 
between shock compression and ramp compression with a high strain rate - but the RH 
equations apply as long as the width of the region of matter where unresolved processes 
occur is constant. Compared with the isentropic states induced by ramp compression in 
a material represented by an EOS, a shock always increases the entropy and hence the 
temperature. With dissipative processes included, the distinction between a ramp and a 
shock may become blurred. 

The RH equations express the conservation of mass, momentum, and energy across a 
moving discontinuity in state. They are usually expressed in terms of the pressure, but are 
readily generalized for materials supporting shear stresses by using the component of stress 
normal to the shock (i.e., parallel with the direction of propagation of the shock), r n : 



where u s is the speed of the shock wave with respect to the material, Au p is the change in 
material speed normal to the shock wave (i.e., parallel to its direction of propagation), and 
subscript refers to the initial state. 

The RH relations can be applied to general material models if a time scale or strain rate 
is imposed, and an orientation chosen for the material with respect to the shock. Shock 
compression in continuum dynamics is almost always uniaxial. 

The RH equations involve only the initial and final states in the material. If a material 




(15) 



(16) 
(17) 
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has properties that depend on the deformation path - such as plastic flow or viscosity 



then physically the detailed shock structure may make a difference [17|. This is a limitation 
of discontinuous shocks in continuum dynamics: it may be addressed as discussed above 
by including dissipative processes and considering ramp compression, if the dissipative pro- 
cesses can be represented adequately in the continuum approximation. Spatially-resolved 
simulations with numerical differentiation to obtain spatial derivatives and forward time 
differencing are usually not capable of representing shock discontinuities directly, and an 
artificial viscosity is used to smear shock compression over a few spatial cells 



The 



trajectory followed by the material in thermodynamic space is a smooth adiabat with dissi- 
pative heating supplied by the artificial viscosity. If plastic work is also included during this 
adiabatic compression, the overall heating for a given compression is greater than from the 
RH equations. To be consistent, plastic flow should be neglected while the artificial viscosity 
is non-zero. This localized disabling of physical processes, particularly time-dependent ones, 
during the passage of the unphysically smeared shock was previously found necessary for 



numerically stable simulations of detonation waves by reactive flow 191 ] . 

Detonation waves are reactive shock waves. Steady planar detonation (the Chapman- 
Jouguet state [^(J) may be calculated using the RH relations, by imposing the condition 
that the material state behind the shock is fully reacted. 

Several numerical methods have been used to solve the RH equations for materials repre- 



sented by an EOS only [21|, [22[. The general RH equations may be solved numerically for a 
given shock compression Ap by varying the specific internal energy e until the normal stress 
from the material model equals that from the RH energy equation, Eq. [TTJ The shock and 
particle speeds are then calculated from Eqs [15] and [161 This numerical method is particu- 
larly convenient for EOS of the form p(p, e), as e may be varied directly. Solutions may still 
be found for general material models using s(e), by which the energy may be varied until 
the solution is found. 

Numerically, the solution was found by bracketing and bisection: 

1. For given compression Ap, take the low-energy end for bracketing as a nearby state 
s_ (e.g. the previous state, of lower compression, on the Hugoniot), compressed adia- 
batically (to state s), and cooled so the specific internal energy is e(s_). 

2. Bracket the desired state: apply successively larger heating increments Ae to s, evolv- 
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ing each trial state internally, until r„(s) from the material model exceeds r n (e — eo) 
from Eq. [T71 

3. Bisect in Ae, evolving each trial state internally, until r n (s) equals r n (e — eo) to the 
desired accuracy. 

As with ramp compression, the independent variable for solution was mass density p, 
and finite steps Ap were taken. Each shock state was calculated independently of the rest, 
so numerical errors did not accumulate along the shock Hugoniot. The accuracy of the 
solution was independent of Ap. A tabular Hugoniot can be calculated by solving over a 
range of p, but again when simulating experimental scenarios it is usually more useful to 
calculate the shock state where one of the other thermodynamic quantities reaches a certain 
value, often that u p and r n match the values from another, simultaneous shock calculation 
for another material - the situation in impact and shock transmission problems, discussed 
below. Specific stop conditions were found by monitoring the quantity of interest until 
bracketed by a finite solution step, then bisecting until the stop condition was satisfied to a 
chosen accuracy. During bisection, each trial calculation was performed as a shock from the 
initial conditions to the trial shock compression. 



C. Accuracy: application to air 

The accuracy of these numerical schemes was tested by comparing with shock and ramp 
compression of a material represented by a perfect gas EOS, 

p=( 7 -l)pe. (18) 

The numerical solution requires a value to be chosen for every parameter in the material 
model, here 7. Air was chosen as an example material, with 7 = 1.4. Air at standard tem- 
perature and pressure has approximately p = 10~ 3 g/cm 3 and e = 0.25MJ/kg. Isentropes 
for the perfect gas EOS have the form 

pp~ 7 = constant, (19) 

and shock Hugoniots have the form 

; ( 7 + l)po-(7-l)p 
11 



The numerical solutions reproduced the principal isentrope and Hugoniot to 10~ 3 % and 0.1% 
respectively, for a compression increment of 1% along the isentrope and a solution tolerance 
of 10~ 6 GPa for each shock state (Fig. [Tj). Over most of the range, the error in the Hugoniot 
was 0.02% or less, only approaching 0.1% near the maximum shock compression. 



IV. COMPLEX BEHAVIOR OF CONDENSED MATTER 

The ability to calculate shock and ramp loci in state space, i.e. as a function of vary- 
ing loading conditions, is particularly convenient for investigating complex aspects of the 
response of condensed matter to dynamic loading. Each locus can be obtained by a single 
series of shock or ramp solutions, rather than having to perform a series of time- and space- 
resolved continuum dynamics simulations, varying the initial or boundary conditions and 
reducing the solution. We consider the calculation of temperature in the scalar EOS, the 
effect of material strength and the effect of phase changes. 



A. Temperature 

The continuum dynamics equations can be closed using a mechanical EOS relating stress 
to mass density, strain, and internal energy. For a scalar EOS, the ideal form to close the 
continuum equations is p(p,e), with s = {p, e} the natural choice for the primitive state 
fields. However, the temperature is needed as a parameter in physical descriptions of many 
contributions to the constitutive response, including plastic flow, phase transitions, and 
chemical reactions. Here, we discuss the calculation of temperature in different forms of the 
scalar EOS. 



1. Density-temperature equations of state 



If the scalar EOS is constructed from its underlying physical contributions for continuum 
dynamics, it may take the form e(p,T), from which p(p,T) can be calculated using the 
second law of thermodynamics 



101. An example is the 'SESAME' form of EOS, based on 

n 

interpolated tabular relations for {p, e}(p, T) [231 ]. A pair of relations {p, e}(p, T) can be 
used as a mechanical EOS by eliminating T, which is equivalent to inverting e(p, T) to find 
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T(p,e), then substituting in p(p, T). For a general e(p, T) relation, for example for the 
SESAME EOS, the inverse can be calculated numerically as required, along an isochore. In 
this way, a {p, e}(p, T) can be used as a p(p, e) EOS. 

Alternatively, the same p(p, T) relation can be used directly with a primitive state field 
including temperature instead of energy: s = {p,T}. The evolution of the state under 
mechanical work then involves the calculation of T(e), i.e. the reciprocal of the specific heat 
capacity, which is a derivative of e(p, T). As this calculation does not require e(p, T) to be 
inverted, it is computationally more efficient to use {p, e}(p, T) EOS with a temperature- 
based, rather than energy-based, state. The main disadvantage is that it is more difficult 
to ensure exact energy conservation as the continuum dynamics equations are integrated in 
time, but any departure from exact conservation is at the level of accuracy of the algorithm 
used to integrate the heat capacity. 

Both structures of EOS have been implemented for material property calculations. Taking 
a SESAME type EOS, thermodynamic loci were calculated with {p, e} or {p, T} primitive 
states, for comparison (Fig. [2]). For a monotonic EOS, the results were indistinguishable 
within differences from forward or reverse interpolation of the tabular relations. When 
the EOS, or the effective surface using a given order of interpolating function, was non- 
monotonic, the results varied greatly because of non-uniqueness when eliminating T for the 
{p, e} primitive state. 

2. Temperature model for mechanical equations of state 

Mechanical EOS are often available as empirical, algebraic relations p(p, e), derived from 
shock data. Temperature can be calculated without altering the mechanical EOS by adding 
a relation T(p, e). While this relation could take any form in principle, one can also follow 
the logic of the Gruneisen EOS, in which the pressure is defined in terms of its deviation 
Ap(p, e — e r ) from a reference curve {p r , e r }(p). Thus temperatures can be calculated by 
reference to a compression curve along which the temperature and specific internal energy 
are known, {T r ,e r }(p), and a specific heat capacity defined as a function of density c„(p). 
In the calculations, this augmented EOS was represented as a 'mechanical-thermal' form 
comprising any p(p, e) EOS plus the reference curves - an example of software inheritance 
and polymorphism. 
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One natural reference curve for temperature is the cold curve, T r = K. The cold curve 
can be estimated from the principal isentrope e(p)| using the estimated density variation 
of the Gruneisen parameter: 

e r (p) = e(p)\ SQ - T c p e a ^M " (21) 



24j . In this work, the principal isentrope was calculated in tabular form from the mechanical 
EOS, using the ramp compression algorithm described above. 

Empirical EOS are calibrated using experimental data. Shock and adiabatic compression 
measurements on strong materials inevitably include elastic-plastic contributions as well as 
the scalar EOS itself. If the elastic-plastic contributions are not taken into account self- 
consistently, the EOS may implicitly include contributions from the strength. A unique 
scalar EOS can be constructed to reproduce the normal stress as a function of compression 
for any unique loading path: shock or adiabat, for a constant or smoothly- varying strain 
rate. Such an EOS would not generally predict the response to other loading histories. The 
EOS and constitutive properties for the materials considered here were constructed self- 
consistently from shock data - this does not mean the models are accurate for other loading 
paths, as neither the EOS nor the strength model includes all the physical terms that real 
materials exhibit. This does not in any case matter for the purposes of demonstrating the 
properties of the numerical schemes. 

This mechanical-thermal procedure was applied to Al using a Gruneisen EOS fitted to the 
same shock data used to calculate the {p, e}(p, T) EOS discussed above [24|. Temperatures 
were in good agreement (Fig. [2]). The mechanical-thermal calculations required a similar 
computational effort to the tabular {p, e}(p, T) EOS with a {p, T} primitive states (and 
were thus much more efficient than the tabular EOS with {p, e} states), and described the 
EOS far more compactly. 

B. Strength 

For dynamic compressions to o(10 GPa) and above, on microsecond time scales, the flow 
stress of solids is often treated as a correction or small perturbation to the scalar EOS. 
Towever, the flow stress has been observed to be much higher on nanosecond time scales 



251 ]. and interactions between elastic and plastic waves may have a significant effect on 
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the compression and wave propagation. The Rankine-Hugoniot equations should be solved 
self-consistently with strength included. 



1. Preferred representation of isotropic strength 

There is an inconsistency in the standard continuum dynamics treatment of scalar (pres- 
sure) and tensor (stress) response. The scalar EOS expresses the pressure p(p, e) as the 
dependent quantity, which is the most convenient form for use in the continuum equations. 
Standard practice is to use sub-Hookean elasticity (hypoelastic form) [16| (Table [IT]), in 
which the state parameters include the stress deviator a, evolved by integration 

a = G(s)e (22) 

where G is the shear modulus and e the strain rate deviator. Thus the isotropic and devia- 
toric contributions to stress are not treated in an equivalent way: the pressure is calculated 
from a local state involving a strain-like parameter (mass density), whereas the stress de- 
viator evolves with the time-derivative of strain. This inconsistency causes problems along 
complicated loading paths because G varies strongly with compression: if a material is sub- 
jected to a shear strain e, then isotropic compression (increasing the shear modulus from 
G to G', leaving e unchanged), then shear unloading to isotropic stress, the true unloading 
strain is — e, whereas the hypoelastic calculation would require a strain of —eG/G'. Using 
Be and the Steinberg- Guinan strength model as an example of the difference between hy- 
poelastic and hyperelastic calculations, consider an initial strain to a flow stress of 0.3 GPa 
followed by isothermal, isotropic compression to 100 GPa,. the strain to unload to a state 
of isotropic stress is 0.20% (hyperelastic) and 0.09% (hypoelastic). The discrepancy arises 
because the hypoelastic model does not increase the deviatoric stress under compression at 
constant deviatoric strain. 

The stress can be considered as a direct response of the material to the instantaneous state 
of elastic strain: a(e, T). This relation can be predicted directly with electronic structure 
calculations of the stress tensor in a solid for a given compression and elastic strain state H] , 
and is a direct generalization of the scalar equation of state. A more consistent representation 
of the state parameters is to use the strain deviator e rather than a, and to calculate a from 
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scratch when required using 

a = G(s)e (23) 

- a hyperelastic formulation. The state parameters are then {p, e, e, e p }. 

The different formulations give different answers when deviatoric strain is accumulated 
at different compressions, in which case the hyperelastic formulation is correct. If the shear 
modulus varies with strain deviator - i.e., for nonlinear elasticity - then the definition of 
G(e) must be adjusted to give the same stress for a given strain. 

Many isotropic strength models use scalar measures of the strain and stress to parame- 
terize work hardening and to apply a yield model of flow stress: 

e=yM- a=y^U\\^\\. (24) 

Inconsistent conventions for equivalent scalar measures have been used by different workers. 
In the present work, the common shock physics convention was used that the flow stress 
component of r n is |y where Y is the flow stress. For consistency with published speeds 
and amplitudes for elastic waves, f e — f a — |, in contrast to other values previously used 



for lower-rate deformation 



26] . In principle, the values of f t and f a do not matter as long as 



the strength parameters were calibrated using the same values then used in any simulations. 



2. Beryllium 

The flow stress measured from laser-driven s 
of micrometers thick is, at around 5-9 GPa 



rock experiments on Be crystals a few tens 



25], much greater than the 0.3-1.3 GPa mea- 



sured on microsecond time scales. A time-dependent crystal plasticity model for Be is being 
developed, and the behavior under dynamic loading depends on the detailed time depen- 
dence of plasticity. Calculations were per formed with the Steinberg-Guinan strength model 



developed for microsecond scale data [24j, and, for the purposes of rough comparison, with 
elastic-perfectly plastic response with a flow stress of 10 GPa. The elastic-perfectly plastic 
model neglected pressure- and work- hardening. 

Calculations were made of the principal adiabat and shock Hugoniot, and of a release 
adiabat from a state on the principal Hugoniot. Calculations were made with and without 
strength. Considering the state trajectories in stress-volume space, it is interesting to note 
that heating from plastic flow may push the adiabat above the Hugoniot, because of the 
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greater heating obtained by integrating along the adiabat compared with jumping from 
the initial to the final state on the Hugoniot (Fig. [3]). Even with an elastic-perfectly plastic 
strength model, the with-strength curves do not lie exactly |y above the strengthless curves, 
because heating from plastic flow contributes an increasing amount of internal energy to the 
EOS as compression increases. 

An important characteristic for the seeding of instabilities by microstructural variations 
in shock response is the shock stress at which an elastic wave does not run ahead of the 
shock. In Be with the high flow stress of nanosecond response, the relation between shock 
and particle speeds is significantly different from the relation for low flow stress (Fig. HJ). For 
low flow stress, the elastic wave travels at 13.2 km/s. A plastic shock travels faster than this 
for pressures greater than 110 GPa, independent of the constitutive model. The speed of a 
plastic shock following the initial elastic wave is similar to the low strength case, because the 
material is already at its flow stress, but the speed of a single plastic shock is appreciably 
higher. 

For compression to a given normal stress, the temperature is significantly higher with 
plastic flow included. The additional heating is particularly striking on the principal adi- 
abat: the temperature departs significantly from the principal isentrope. Thus ramp-wave 
compression of strong materials may lead to significant levels of heating, contrary to com- 



mon assumptions of small temperature increases [27]]. Plastic flow is largely irreversible, so 
heating occurs on unloading as well as loading. Thus, on adiabatic release from a shock- 
compressed state, additional heating occurs compared with the no-strength case. These 
levels of heating are important as shock or release melting may occur at a significantly lower 
shock pressure than would be expected ignoring the effect of strength. (Fig. 0) 



C. Phase changes 

An important property of condensed matter is phase changes, including solid-solid poly- 
morphism and solid-liquid. An equilibrium phase diagram can be represented as a single 
overall EOS surface as before. Multiple, competing phases with kinetics for each phase trans- 
formation can be represented conveniently using the structure described above for general 
material properties, for example by describing the local state as a set of volume fractions 
/, of each possible simple-EOS phase, with transition rates and equilibration among them. 
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This model is described in more detail elsewhere 19J. However, it is interesting to investi- 
gate the robustness of the numerical scheme for calculating shock Hugoniots when the EOS 
has the discontinuities in value and gradient associated with phase changes. 

The EOS of molten metal, and the solid-liquid phase transition, can be represented to a 
reasonable approximation as an adjustment to the EOS of the solid: 



where 



Ptwo-phase(P> e ) =Psolid(^ g ) 



(25) 



(26) 



e : T(p,e) < T m (p) 

e - Ah m : Ah m = c v (p, e) [T(p, e) - T m (p)} < Ah m 
e — Ah m : otherwise 
and Ah m is the specific latent heat of fusion. Taking the EOS and a modified Lindemann 



melting curve for Al 24| . and using Ah m = 0.397 MJ/kg, the shock Hugoniot algorithm was 



found to operate stably across the phase transition (Fig. [6]). 



V. COMPOSITE LOADING PATHS 



Given methods to calculate shock and adiabatic loading paths from arbitrary initial 
states, a considerable variety of experimental scenarios can be treated from the interaction 
of loading or unloading waves with interfaces between different materials, in planar geometry 
for uniaxial compression. The key physical constraint is that, if two dissimilar materials are 
to remain in contact after an interaction such as an impact or the passage of a shock, the 
normal stress r n and particle speed u p in both materials must be equal on either side of the 
interface. The change in particle speed and stress normal to the waves were calculated above 
for compression waves running in the direction of increasing spatial ordinate (left to right). 
Across an interface, the sense is reversed for the material at the left. Thus a projectile 
impacting a stationary target to the right is decelerated from its initial speed by the shock 
induced by impact. 

The general problem at an interface can be analyzed by considering the states at the 
instant of first contact - on impact, or when a shock traveling through a sandwich of ma- 
terials first reaches the interface. The initial states are {ui, sf, u r , s r }. The final states are 
{uj, s'f, Uj, r' r } where Uj is the joint particle speed, r n (sj) = r n (s' r ), and s[ is connected to Sj 
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by either a shock or an adiabat, starting at the appropriate initial velocity and stress, and 
with orientation given by the side of the system each material occurs on. Each type of wave 
is considered in turn, looking for an intersection in the u p — r n plane. Examples of these 
wave interactions are the impact of a projectile with a stationary target (Fig. [7]), release of a 
shock state at a free surface or a material (e.g. a window) of lower shock impedance (hence 
reflecting a release wave into the shocked material - Fig. [8]), reshocking at a surface with a 
material of higher shock impedance (Fig. [8]), or tension induced as materials try to separate 
in opposite directions when joined by a bonded interface (Fig. [9]). Each of these scenarios 
may occur in turn following the impact of a projectile with a target: if the target is layered 
then a shock is transmitted across each interface with a release or a reshock reflected back, 
depending on the materials; release ultimately occurs at the rear of the projectile and the 
far end of the target, and the oppositely-moving release waves subject the projectile and 
target to tensile stresses when they interact (Fig. [TO]) . 

As an illustration of combining shock and ramp loading calculations, consider the problem 
of an Al projectile, initially traveling at 3.6 km/s, im pac ting a stationary, composite target 



comprising a Mo sample and a LiF release window [28j, [29] . The shock and release states were 
calculated using published material properties 24] . The initial shock state was calculated to 
have a normal stress of 63.9 GPa. On reaching the LiF, the shock was calculated to transmit 
at 27.1 GPa, reflecting as a release in the Mo. These stresses match the continuum dynamics 
simulation to within 0.1 GPa in the Mo and 0.3 GPa in the LiF, using the same material 
properties (Fig. [TT1) . The associated wave and particle speeds match to a similar accuracy; 
wave speeds are much more difficult to extract from the continuum dynamics simulation. 
An extension of this analysis can be used to calculate the interaction of oblique shocks 

r 

with an interface I3( 



VI. CONCLUSIONS 



A general formulation was developed to represent material models for applications in 
dynamic loading, suitable for software implementation in object-oriented programming lan- 
guages. Numerical methods were devised to calculate the response of matter represented 
by the general material models to shock and ramp compression, and ramp decompression, 
by direct evaluation of the thermodynamic pathways for these compressions rather than 
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spatially-resolved simulations. This approach is a generalization of earlier work on solutions 
for materials represented by a scalar equation of state. The numerical methods were found 
to be flexible and robust: capable of application to materials with very different properties. 
The numerical solutions matched analytic results to a high accuracy. 

Care was needed with the interpretation of some types of physical response, such as plas- 
tic flow, when applied to deformation at high strain rates. The underlying time-dependence 
of processes occurring during deformation should be taken into account. The actual history 
of loading and heating experienced by material during the passage of a shock may influence 
the final state - this history is not captured in the continuum approximation to material 
dynamics, where shocks are treated as discontinuities. Thus care is also needed in spa- 
tially resolved simulations when shocks are modeled using artificial viscosity to smear them 
unphysically over a finite thickness. 

Calculations were shown to demonstrate the operation of the algorithms for shock and 
ramp compression with material models representative of complex solids including strength 
and phase transformations. 

The basic ramp and shock solution methods were coupled to solve for composite defor- 
mation paths, such as shock-induced impacts, and shock interactions with a planar interface 
between different materials. Such calculations capture much of the physics of typical ma- 
terial dynamics experiments, without requiring spatially-resolving simulations. The results 
of direct solution of the relevant shock and ramp loading conditions were compared with 
hydrocode simulations, showing complete consistency. 
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TABLE I: Interface to material models required for explicit forward-time continuum dynamics 
simulations. 



purpose 


interface calls 


program set-up 


read/write material data 


continuum dynamics equations 


stress (state) 


time step control 


sound speed (state) 


evolution of state (deformation) c/(state)/dt(state,grad'iZ) 


evolution of state (heating) 


d(state) / dt ( state, e) 


internal evolution of state 


d(staie)/dt 


manipulation of states 


create and delete 




add states 




multiply state by a scalar 




check for self-consistency 



Parentheses in the interface calls denote functions, e.g. "stress (state)" for "stress as a function of 
the instantaneous, local state." The evolution functions are shown in the operator-split structure 
that is most robust for explicit, forward-time numerical solutions and can also be used for 
calculations of the shock Hugoniot and ramp compression. Checks for self-consistency include 
that mass density is positive, volume or mass fractions of components of a mixture add up to one, 

etc. 
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TABLE II: Examples of types of material model, distinguished by different structures in the state 
vector. 



model state vector effect of mechanical strain 

s s m (s,giadu) 

mechanical equation of state p, e —pdivu, — pdivu/ p 

thermal equation of state p, T —pdivu, —pdivu/ pc v 

heterogeneous mixture {p,e,f v }i {—pdivu, —pdivu/ p,0}i 

homogeneous mixture p,T,{f m }i {—pdivu, —pdivu/ pc v ,0i 

traditional deviatoric strength p, e, a, e p —pdivu, P^ v "+/pll CTe pll ^ Ge e , 



The symbols are p: mass density; e: specific internal energy, T: temperature, /„: volume fraction, 
f m : mass fraction, a: stress deviator, f p : fraction of plastic work converted to heat, grad« p : 
plastic part of velocity gradient, G: shear modulus, e eiP : elastic and plastic parts of strain rate 
deviator, e p : scalar equivalent plastic strain, f e : factor in effective strain magnitude. Reacting 
solid explosives can be represented as heterogeneous mixtures, one component being the reacted 
products; reaction, a process of internal evolution, transfers material from unreacted to reacted 
components. Gas-phase reaction can be represented as a homogeneous mixture, reactions 
transferring mass between components representing different types of molecule. Symmetric 
tensors such as the stress deviator are represented more compactly by their 6 unique upper 
triangular components, e.g. using Voigt notation. 
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TABLE III: Outline hierarchy of material models, illustrating the use of polymorphism (in the 
object-oriented programming sense). 



material (or state) type model type 



thermal equation of state 



mechanical equation of state polytropic, Griineisen, energy-based 

Jones- Wilkins-Lee, (p,T) table, etc 
temperature-based Jones- Wilkins- 
Lee, quasiharmonic, (p, T) table, 
etc 

modified polytropic, reactive Jones- 
Wilkins-Lee 
Cochran-Banner 
elastic-plastic, 
Steinberg-Lund. 
Wallace, etc 

mixing and reaction models 
equilibration and reaction models 



reactive equation of state 
spall 

deviatoric stress 



S teinb erg- Guinan , 
Preston- Tonks- 



homogeneous mixture 
heterogeneous mixture 



Continuum dynamics programs can refer to material properties as an abstract 'material type' 
with an abstract material state. The actual type of a material (e.g. mechanical equation of 
state), the specific model type (e.g. polytropic), and the state of material of that type are all 
handled transparently by the object-oriented software structure. 
The reactive equation of state has an additional state parameter A, and the software operations 

are defined by extending those of the mechanical equation of state. Spalling materials can be 
represented by a solid state plus a void fraction f v , with operations defined by extending those of 
the solid material. Homogeneous mixtures are defined as a set of thermal equations of state, and 
the state is the set of states and mass fractions for each. Heterogeneous mixtures are defined as a 
set of 'pure' material properties of any type, and the state is the set of states for each component 

plus its volume fraction. 
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FIG. 1: Principal isentrope and shock Hugoniot for air (perfect gas): numerical calculations for 
general material models, compared with analytic solutions. 
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temperature (K) 



FIG. 2: Shock Hugoniot for Al in pressure-temperature space, for different representations of the 
equation of state. 
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FIG. 3: Principal adiabat and shock Hugoniot for Be in normal stress-compression space, neglecting 
strength (dashed), for Steinberg-Guinan strength (solid), and for elastic-perfectly plastic with 
Y = lOGPa (dotted). 
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FIG. 4: Principal adiabat and shock Hugoniot for Be in shock speed-normal stress space, neglecting 
strength (dashed), for Steinberg-Guinan strength (solid), and for elastic-perfectly plastic with 
Y = 10 GPa (dotted). 
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FIG. 5: Principal adiabat, shock Hugoniot, and release adiabat for Be in normal stress-temperature 
space, neglecting strength (dashed), for Steinberg-Guinan strength (solid), and for elastic-perfectly 
plastic with Y = 10 GPa (dotted). 
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FIG. 6: Demonstration of shock Hugoniot solution across a phase boundary: shock-melting of Al, 
for different initial porosities. 
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FIG. 7: Wave interactions for the impact of a flat projectile moving from left to right with 
stationary target. Dashed arrows are a guide to the sequence of states. For a projectile movin 
from right to left, the construction is the mirror image reflected in the normal stress axis. 
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target release at free surface 



FIG. 8: Wave interactions for the release of a shocked state (shock moving from left to right) into 
a stationary 'window' material to its right. The release state depends whether the window has 
a higher or lower shock impedance than the shocked material. Dashed arrows are a guide to the 
sequence of states. For a shock moving from right to left, the construction is the mirror image 
reflected in the normal stress axis. 
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FIG. 9: Wave interactions for the release of a shocked state by tension induced as materials try- 
to separate in opposite directions when joined by a bonded interface. Material damage, spall, and 
separation are neglected: the construction shows the maximum tensile stress possible. For general 
material properties, e.g. if plastic flow is included, the state of maximum tensile stress is not just 
the negative of the initial shock state. Dashed arrows are a guide to the sequence of states. The 
graph shows the initial state after an impact by a projectile moving from right to left; for a shock 
moving from right to left, the construction is the mirror image reflected in the normal stress axis. 
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FIG. 10: Schematic of uniaxial wave interactions induced by the impact of a flat projectile with a 
composite target. 
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FIG. 11: Hydrocode simulation of Al projectile at 3.6km/s impacting a Mo target with a LiF 
release window, 1.1 /is after impact. Structures on the waves are elastic precursors. 
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